SRCD 2019 Symposium: Religious & metaphysical concepts (Srinivasan)

# load required libraries
library(tidyverse)
library(langcog) # source: https://github.com/langcog/langcog-package
library(psych)
library(lme4)
# set theme for ggplots
theme_set(theme_bw())

Data preparation

d_all <- d_us_ad_pilot %>% rownames_to_column("subid") %>%
  full_join(d_gh_ch %>% rownames_to_column("subid")) %>%
  full_join(d_th_ad %>% rownames_to_column("subid")) %>%
  full_join(d_th_ch %>% rownames_to_column("subid")) %>%
  full_join(d_vt_ad %>% rownames_to_column("subid")) %>%
  full_join(d_vt_ch %>% rownames_to_column("subid")) %>%
  column_to_rownames("subid")

Joining, by = c("subid", "get angry", "choose what to do", "figure out how to do things", "feel guilty", "feel happy", "hear things", "get hungry", "get hurt feelings", "feel love", "feel pain", "pray", "feel proud", "remember things", "feel sad", "feel scared", "sense when things are far away", "sense temperatures", "feel shy", "feel sick, like when you feel like you might vomit", "smell things", "think about things", "feel tired")
Joining, by = c("subid", "get angry", "choose what to do", "figure out how to do things", "feel guilty", "feel happy", "hear things", "get hungry", "get hurt feelings", "feel love", "feel pain", "pray", "feel proud", "remember things", "feel sad", "feel scared", "sense when things are far away", "sense temperatures", "feel shy", "feel sick, like when you feel like you might vomit", "smell things", "think about things", "feel tired")
Joining, by = c("subid", "get angry", "choose what to do", "figure out how to do things", "feel guilty", "feel happy", "hear things", "get hungry", "get hurt feelings", "feel love", "feel pain", "pray", "feel proud", "remember things", "feel sad", "feel scared", "sense when things are far away", "sense temperatures", "feel shy", "feel sick, like when you feel like you might vomit", "smell things", "think about things", "feel tired")
Joining, by = c("subid", "get angry", "choose what to do", "figure out how to do things", "feel guilty", "feel happy", "hear things", "get hungry", "get hurt feelings", "feel love", "feel pain", "pray", "feel proud", "remember things", "feel sad", "feel scared", "sense when things are far away", "sense temperatures", "feel shy", "feel sick, like when you feel like you might vomit", "smell things", "think about things", "feel tired")
Joining, by = c("subid", "get angry", "choose what to do", "figure out how to do things", "feel guilty", "feel happy", "hear things", "get hungry", "get hurt feelings", "feel love", "feel pain", "pray", "feel proud", "remember things", "feel sad", "feel scared", "sense when things are far away", "sense temperatures", "feel shy", "feel sick, like when you feel like you might vomit", "smell things", "think about things", "feel tired")

Shared conceptual structure

Pooling all participants from all sites together into a common factor structure.

Parallel analysis

How many factors to retain?

reten_all_PA <- fa.parallel(d_all, plot = F); reten_all_PA

Parallel analysis suggests that the number of factors =  4  and the number of components =  3 
Call: fa.parallel(x = d_all, plot = F)
Parallel analysis suggests that the number of factors =  4  and the number of components =  3 

 Eigen Values of 
  Original factors Simulated data Original components simulated data
1             8.53           0.33                9.13           1.29
2             1.28           0.26                1.92           1.24
3             0.62           0.22                1.26           1.21
4             0.20           0.19                0.83           1.17

reten_all_par <- reten_all_PA$nfact

What are these factors?

efa_all_par <- fa(d_all, nfactors = reten_all_par, rotate = "oblimin",
                  scores = "tenBerge", impute = "median")

convergence not obtained in GPFoblq. 1000 iterations used.

heatmap_fun(efa_all_par) + 
  labs(title = "Parallel Analysis")

Joining, by = "capacity"
Joining, by = "factor"

Which capacities are attributed to which targets?

scoresplot_fun(efa_all_par, target = "all") + 
  labs(title = "Parallel Analysis")

Ignoring unknown aesthetics: y

scoresplot_fun(efa_all_par, target = c("ghosts", "God", "children")) + 
  labs(title = "Parallel Analysis")

the condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedIgnoring unknown aesthetics: y

itemsplot_fun(efa_all_par, target = c("ghosts", "God", "children")) + 
  labs(title = "Parallel Analysis")

the condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedJoining, by = "capacity"

|============================                                 | 47% ~2 s remaining     
|=============================                                | 48% ~2 s remaining     
|==============================                               | 51% ~2 s remaining     
|===============================                              | 52% ~2 s remaining     
|================================                             | 54% ~2 s remaining     
|=================================                            | 55% ~2 s remaining     
|===================================                          | 58% ~2 s remaining     
|====================================                         | 59% ~2 s remaining     
|=====================================                        | 61% ~2 s remaining     
|======================================                       | 63% ~1 s remaining     
|=======================================                      | 65% ~1 s remaining     
|========================================                     | 66% ~1 s remaining     
|=========================================                    | 68% ~1 s remaining     
|==========================================                   | 69% ~1 s remaining     
|===========================================                  | 71% ~1 s remaining     
|============================================                 | 72% ~1 s remaining     
|=============================================                | 74% ~1 s remaining     
|==============================================               | 76% ~1 s remaining     
|===============================================              | 78% ~1 s remaining     
|================================================             | 80% ~1 s remaining     
|=================================================            | 81% ~1 s remaining     
|==================================================           | 83% ~1 s remaining     
|===================================================          | 84% ~1 s remaining     
|===================================================          | 85% ~1 s remaining     
|====================================================         | 86% ~1 s remaining     
|=====================================================        | 87% ~0 s remaining     
|=====================================================        | 88% ~0 s remaining     
|=====================================================        | 88% ~0 s remaining     
|======================================================       | 89% ~0 s remaining     
|======================================================       | 89% ~0 s remaining     
|======================================================       | 90% ~0 s remaining     
|=======================================================      | 92% ~0 s remaining     
|========================================================     | 93% ~0 s remaining     
|========================================================     | 93% ~0 s remaining     
|=========================================================    | 94% ~0 s remaining     
|==========================================================   | 95% ~0 s remaining     
|==========================================================   | 96% ~0 s remaining     
|===========================================================  | 98% ~0 s remaining     
|============================================================ | 99% ~0 s remaining     
|=============================================================|100% ~0 s remaining     

Joining, by = c("factor", "capacity", "order")

Minimizing BIC

How many factors to retain?

reten_all_BIC <- VSS(d_all, plot = F); reten_all_BIC

Very Simple Structure
Call: vss(x = x, n = n, rotate = rotate, diagonal = diagonal, fm = fm, 
    n.obs = n.obs, plot = plot, title = title, use = use, cor = cor)
VSS complexity 1 achieves a maximimum of 0.88  with  1  factors
VSS complexity 2 achieves a maximimum of 0.92  with  2  factors

The Velicer MAP achieves a minimum of 0.01  with  3  factors 
BIC achieves a minimum of  -645.05  with  3  factors
Sample Size adjusted BIC achieves a minimum of  -196.92  with  5  factors

Statistics by number of factors 
  vss1 vss2   map dof chisq     prob sqresid  fit RMSEA  BIC SABIC complex eChisq
1 0.88 0.00 0.018 209  1964 1.7e-282    11.0 0.88 0.096  539  1203     1.0   2336
2 0.56 0.92 0.011 188   899  4.0e-93     7.5 0.92 0.065 -383   214     1.5    696
3 0.43 0.82 0.011 168   501  9.0e-35     6.2 0.93 0.047 -645  -112     1.8    265
4 0.43 0.81 0.014 149   383  1.5e-22     5.7 0.94 0.042 -633  -160     1.9    184
5 0.40 0.73 0.018 131   280  6.9e-13     5.4 0.94 0.036 -613  -197     2.1    134
6 0.40 0.73 0.023 114   227  1.6e-09     5.1 0.95 0.033 -550  -188     2.2    103
7 0.40 0.68 0.027  98   182  5.7e-07     4.8 0.95 0.031 -486  -175     2.3     75
8 0.37 0.60 0.033  83   134  3.3e-04     4.5 0.95 0.026 -432  -168     2.5     52
   SRMR eCRMS eBIC
1 0.074 0.078  911
2 0.041 0.045 -586
3 0.025 0.029 -880
4 0.021 0.026 -832
5 0.018 0.024 -759
6 0.016 0.022 -674
7 0.013 0.020 -594
8 0.011 0.019 -514

reten_all_bic <- data.frame(reten_all_BIC$vss.stats %>% rownames_to_column("nfact") %>% top_n(-1, BIC))$nfact %>% as.numeric()

What are these factors?

efa_all_bic <- fa(d_all, nfactors = reten_all_bic, rotate = "oblimin",
                  scores = "tenBerge", impute = "median")

heatmap_fun(efa_all_bic) + 
  labs(title = "Minimizing BIC")

Joining, by = "capacity"
Joining, by = "factor"

  # labs(title = "Figure 1: Shared conceptual structure")

Which capacities are attributed to which targets?

scoresplot_fun(efa_all_bic, target = "all") + #, highlight = "supernatural") + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1,
                                   color = c(rep("black", 8),
                                             rep("#984ea3", 2)),
                                   face = c(rep("plain", 8),
                                            rep("bold", 2)))) +
  scale_fill_manual(values = c("#e41a1c", "#4daf4a", "#377eb8")) +
  labs(title = "Minimizing BIC")

Ignoring unknown aesthetics: yScale for 'fill' is already present. Adding another scale for 'fill', which will
replace the existing scale.

  # labs(title = "Figure 2: Factor scores")

scoresplot_fun(efa_all_bic, target = c("ghosts", "God", "children")) + 
  scale_fill_manual(values = c("#e41a1c", "#4daf4a", "#377eb8")) +
  labs(title = "Minimizing BIC")

the condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedIgnoring unknown aesthetics: yScale for 'fill' is already present. Adding another scale for 'fill', which will
replace the existing scale.

itemsplot_fun(efa_all_bic, target = c("ghosts", "God", "children")) + 
  labs(title = "Minimizing BIC")

the condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedJoining, by = "capacity"

|===============================================              | 78% ~1 s remaining     
|=================================================            | 81% ~0 s remaining     
|==================================================           | 83% ~0 s remaining     
|===================================================          | 85% ~0 s remaining     
|=====================================================        | 87% ~0 s remaining     
|======================================================       | 89% ~0 s remaining     
|=======================================================      | 91% ~0 s remaining     
|========================================================     | 93% ~0 s remaining     
|========================================================     | 93% ~0 s remaining     
|=========================================================    | 95% ~0 s remaining     
|==========================================================   | 95% ~0 s remaining     
|==========================================================   | 96% ~0 s remaining     
|===========================================================  | 98% ~0 s remaining     
|============================================================ | 99% ~0 s remaining     

Joining, by = c("factor", "capacity", "order")

Preset criteria

How many factors to retain?

reten_all_k <- reten_fun(d_all, rot_type = "oblimin"); reten_all_k

[1] 2

What are these factors?

efa_all_k <- fa(d_all, nfactors = reten_all_k, rotate = "oblimin",
                  scores = "tenBerge", impute = "median")

heatmap_fun(efa_all_k) + 
  labs(title = "Preset criteria (Weisman et al., 2017)")

Joining, by = "capacity"
Joining, by = "factor"

Which capacities are attributed to which targets?

scoresplot_fun(efa_all_k, target = "all") + 
  labs(title = "Preset criteria (Weisman et al., 2017)")

Ignoring unknown aesthetics: y

scoresplot_fun(efa_all_k, target = c("ghosts", "God", "children")) + 
  labs(title = "Preset criteria (Weisman et al., 2017)")

the condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedIgnoring unknown aesthetics: y

itemsplot_fun(efa_all_k, target = c("ghosts", "God", "children")) + 
  labs(title = "Preset criteria (Weisman et al., 2017)")

the condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedJoining, by = "capacity"

|==============================================               | 76% ~1 s remaining     
|===============================================              | 77% ~1 s remaining     
|===============================================              | 78% ~1 s remaining     
|================================================             | 80% ~1 s remaining     
|==================================================           | 83% ~0 s remaining     
|===================================================          | 85% ~0 s remaining     
|====================================================         | 87% ~0 s remaining     
|=====================================================        | 88% ~0 s remaining     
|=======================================================      | 90% ~0 s remaining     
|========================================================     | 93% ~0 s remaining     
|=========================================================    | 95% ~0 s remaining     
|==========================================================   | 95% ~0 s remaining     
|===========================================================  | 98% ~0 s remaining     
|============================================================ | 99% ~0 s remaining     

Joining, by = c("factor", "capacity", "order")

Ghana: children

Parallel analysis

How many factors to retain?

reten_gh_ch_PA <- fa.parallel(d_gh_ch, plot = F); reten_gh_ch_PA

Parallel analysis suggests that the number of factors =  2  and the number of components =  2 
Call: fa.parallel(x = d_gh_ch, plot = F)
Parallel analysis suggests that the number of factors =  2  and the number of components =  2 

 Eigen Values of 
  Original factors Simulated data Original components simulated data
1             7.86           0.88                8.48           1.78
2             1.96           0.68                2.62           1.65

reten_gh_ch_par <- reten_gh_ch_PA$nfact

What are these factors?

efa_gh_ch_par <- fa(d_gh_ch, nfactors = reten_gh_ch_par, rotate = "oblimin",
                  scores = "tenBerge", impute = "median")

heatmap_fun(efa_gh_ch_par) + 
  labs(title = "Parallel Analysis")

Joining, by = "capacity"
Joining, by = "factor"

Which capacities are attributed to which targets?

scoresplot_fun(efa_gh_ch_par, target = "all") + 
  labs(title = "Parallel Analysis")

Ignoring unknown aesthetics: y

scoresplot_fun(efa_gh_ch_par, target = c("ghosts", "God", "children")) + 
  labs(title = "Parallel Analysis")

the condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedIgnoring unknown aesthetics: y

itemsplot_fun(efa_gh_ch_par, target = c("ghosts", "God", "children")) + 
  labs(title = "Parallel Analysis")

the condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedJoining, by = "capacity"

|============================================                 | 73% ~1 s remaining     
|=============================================                | 75% ~1 s remaining     
|==============================================               | 77% ~1 s remaining     
|===============================================              | 79% ~1 s remaining     
|================================================             | 80% ~1 s remaining     
|=================================================            | 81% ~1 s remaining     
|==================================================           | 83% ~0 s remaining     
|====================================================         | 86% ~0 s remaining     
|=====================================================        | 88% ~0 s remaining     
|======================================================       | 89% ~0 s remaining     
|=======================================================      | 91% ~0 s remaining     
|========================================================     | 93% ~0 s remaining     
|==========================================================   | 95% ~0 s remaining     
|===========================================================  | 97% ~0 s remaining     
|===========================================================  | 98% ~0 s remaining     
|=============================================================|100% ~0 s remaining     

Joining, by = c("factor", "capacity", "order")

Minimizing BIC

How many factors to retain?

reten_gh_ch_BIC <- VSS(d_gh_ch, plot = F); reten_gh_ch_BIC

Very Simple Structure
Call: vss(x = x, n = n, rotate = rotate, diagonal = diagonal, fm = fm, 
    n.obs = n.obs, plot = plot, title = title, use = use, cor = cor)
VSS complexity 1 achieves a maximimum of 0.83  with  1  factors
VSS complexity 2 achieves a maximimum of 0.91  with  2  factors

The Velicer MAP achieves a minimum of 0.02  with  2  factors 
BIC achieves a minimum of  -597.63  with  2  factors
Sample Size adjusted BIC achieves a minimum of  -57.2  with  8  factors

Statistics by number of factors 
  vss1 vss2   map dof chisq    prob sqresid  fit RMSEA  BIC SABIC complex eChisq  SRMR
1 0.83 0.00 0.034 209   645 5.4e-46    14.8 0.83 0.124 -401 260.6     1.0    824 0.109
2 0.62 0.91 0.017 188   343 3.7e-11     8.1 0.91 0.080 -598  -2.7     1.4    199 0.054
3 0.58 0.89 0.019 168   282 9.2e-08     7.0 0.92 0.074 -559 -27.4     1.7    145 0.046
4 0.51 0.85 0.022 149   247 7.5e-07     6.1 0.93 0.073 -498 -26.8     1.8    114 0.041
5 0.48 0.81 0.025 131   208 2.3e-05     5.4 0.94 0.070 -448 -33.2     2.0     85 0.035
6 0.46 0.80 0.030 114   163 1.9e-03     4.9 0.94 0.061 -408 -47.1     2.1     64 0.031
7 0.43 0.76 0.034  98   132 1.3e-02     4.3 0.95 0.057 -358 -48.3     2.2     48 0.027
8 0.34 0.65 0.040  83    95 1.7e-01     3.9 0.96 0.042 -320 -57.2     2.5     34 0.022
  eCRMS eBIC
1 0.115 -222
2 0.060 -742
3 0.054 -696
4 0.051 -632
5 0.047 -570
6 0.043 -506
7 0.041 -442
8 0.037 -381

reten_gh_ch_bic <- data.frame(reten_gh_ch_BIC$vss.stats %>% rownames_to_column("nfact") %>% top_n(-1, BIC))$nfact %>% as.numeric()

What are these factors?

efa_gh_ch_bic <- fa(d_gh_ch, nfactors = reten_gh_ch_bic, rotate = "oblimin",
                  scores = "tenBerge", impute = "median")

heatmap_fun(efa_gh_ch_bic) + 
  labs(title = "Minimizing BIC")

Joining, by = "capacity"
Joining, by = "factor"

  # labs(title = "Figure 1: Shared conceptual structure")

Which capacities are attributed to which targets?

scoresplot_fun(efa_gh_ch_bic, target = "all") + #, highlight = "supernatural") + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1,
                                   color = c(rep("black", 8),
                                             rep("#984ea3", 2)),
                                   face = c(rep("plain", 8),
                                            rep("bold", 2)))) +
  scale_fill_manual(values = c("#e41a1c", "#4daf4a", "#377eb8")) +
  labs(title = "Minimizing BIC")

Ignoring unknown aesthetics: yScale for 'fill' is already present. Adding another scale for 'fill', which will
replace the existing scale.

  # labs(title = "Figure 2: Factor scores")

scoresplot_fun(efa_gh_ch_bic, target = c("ghosts", "God", "children")) + 
  scale_fill_manual(values = c("#e41a1c", "#4daf4a", "#377eb8")) +
  labs(title = "Minimizing BIC")

the condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedIgnoring unknown aesthetics: yScale for 'fill' is already present. Adding another scale for 'fill', which will
replace the existing scale.

itemsplot_fun(efa_gh_ch_bic, target = c("ghosts", "God", "children")) + 
  labs(title = "Minimizing BIC")

the condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedJoining, by = "capacity"

|====================================================         | 86% ~0 s remaining     
|=====================================================        | 87% ~0 s remaining     
|======================================================       | 89% ~0 s remaining     
|=======================================================      | 92% ~0 s remaining     
|=========================================================    | 94% ~0 s remaining     
|==========================================================   | 96% ~0 s remaining     
|============================================================ | 98% ~0 s remaining     
|=============================================================|100% ~0 s remaining     

Joining, by = c("factor", "capacity", "order")

Preset criteria

How many factors to retain?

reten_gh_ch_k <- reten_fun(d_gh_ch, rot_type = "oblimin"); reten_gh_ch_k

[1] 2

What are these factors?

efa_gh_ch_k <- fa(d_gh_ch, nfactors = reten_gh_ch_k, rotate = "oblimin",
                  scores = "tenBerge", impute = "median")

heatmap_fun(efa_gh_ch_k) + 
  labs(title = "Preset criteria (Weisman et al., 2017)")

Joining, by = "capacity"
Joining, by = "factor"

Which capacities are attributed to which targets?

scoresplot_fun(efa_gh_ch_k, target = "all") + 
  labs(title = "Preset criteria (Weisman et al., 2017)")

Ignoring unknown aesthetics: y

scoresplot_fun(efa_gh_ch_k, target = c("ghosts", "God", "children")) + 
  labs(title = "Preset criteria (Weisman et al., 2017)")

the condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedIgnoring unknown aesthetics: y

itemsplot_fun(efa_gh_ch_k, target = c("ghosts", "God", "children")) + 
  labs(title = "Preset criteria (Weisman et al., 2017)")

the condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedthe condition has length > 1 and only the first element will be usedJoining, by = "capacity"

|===============================================              | 78% ~1 s remaining     
|================================================             | 79% ~1 s remaining     
|=================================================            | 81% ~1 s remaining     
|==================================================           | 83% ~0 s remaining     
|===================================================          | 84% ~0 s remaining     
|====================================================         | 86% ~0 s remaining     
|=====================================================        | 88% ~0 s remaining     
|======================================================       | 89% ~0 s remaining     
|=======================================================      | 91% ~0 s remaining     
|========================================================     | 93% ~0 s remaining     
|=========================================================    | 94% ~0 s remaining     
|==========================================================   | 96% ~0 s remaining     
|===========================================================  | 98% ~0 s remaining     
|=============================================================|100% ~0 s remaining     

Joining, by = c("factor", "capacity", "order")